This tutorial extends the 1-pool Bayesian tutorial to a 3-pool series model. If you haven’t gone through the 1-pool tutorial yet, start there – it introduces Bayesian inference, likelihood functions, priors, and MCMC from scratch.
In the 1-pool tutorial, we learned that:
k
and carbon input) control the modelNow we move to 6 parameters and 6 data points. The problem is formally just as determined, but the parameter space is much harder to explore and the parameters can trade off against each other in complex ways.
library(SoilR)
library(BayesianTools)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
set.seed(42)
load("soilRmcmc.RData")
We reuse the same run_and_summarize helper from the
1-pool tutorial, extended to handle any number of parameters.
# 'prior' can be any BayesianTools prior object (uniform, truncated normal, custom).
# If NULL, a uniform prior is created from prior_lower/prior_upper.
run_and_summarize <- function(likelihood_fn, prior = NULL,
prior_lower = NULL, prior_upper = NULL,
par_names,
iterations = 20000,
burn_in = 5000, thin = 10) {
if (is.null(prior)) {
prior <- createUniformPrior(lower = prior_lower, upper = prior_upper)
}
bt_setup <- createBayesianSetup(
likelihood = likelihood_fn,
prior = prior,
names = par_names
)
mcmc_out <- runMCMC(
bayesianSetup = bt_setup,
sampler = "DREAMzs",
settings = list(iterations = iterations, message = FALSE)
)
raw_chains <- getSample(mcmc_out, start = 1, coda = TRUE, parametersOnly = FALSE)
chain_df <- do.call(rbind, lapply(seq_along(raw_chains), function(i) {
ch <- as.data.frame(as.matrix(raw_chains[[i]]))
n_par <- length(par_names)
colnames(ch)[1:n_par] <- par_names
ch$log_posterior <- ch[, n_par + 1]
ch <- ch[, c(par_names, "log_posterior")]
ch$iteration <- seq_len(nrow(ch))
ch$chain <- factor(i)
ch
}))
samples <- getSample(mcmc_out, start = burn_in, thin = thin)
colnames(samples) <- par_names
posterior_df <- as.data.frame(samples)
map_est <- MAP(mcmc_out, start = burn_in)$parametersMAP
names(map_est) <- par_names
summary_tbl <- posterior_df %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
group_by(parameter) %>%
summarise(
median = median(value),
lower_90 = quantile(value, 0.05),
upper_90 = quantile(value, 0.95),
.groups = "drop"
)
list(
mcmc_out = mcmc_out,
chain_df = chain_df,
posterior_df = posterior_df,
map = map_est,
summary = summary_tbl,
burn_in = burn_in,
thin = thin
)
}
The CiPEHR soil profile has three distinct layers: surface organic, deep organic, and mineral. These layers have very different properties – the surface is fresh, carbon-rich litter that decomposes quickly, while the mineral layer is old, stabilized carbon that turns over on centennial timescales.
When we aggregated everything into a single pool in the 1-pool tutorial, we lost this structure. The bulk \(\Delta^{14}\)C is a weighted average that hides how the bomb spike propagates differently through each layer.
Let’s look at the per-layer data.
dat_layers <- datComb %>%
filter(year == 2009 | treatment == "Control") %>%
mutate(layer_label = factor(layer,
levels = c("so", "do", "min"),
labels = c("Surface organic", "Deep organic", "Mineral")))
knitr::kable(
dat_layers %>%
select(year, layer_label, C_stock, C_stock_sd, C14, C_14_sd) %>%
arrange(year, layer_label),
digits = 1,
col.names = c("Year", "Layer", "C stock (kg/m2)", "C stock SD",
"Delta14C", "Delta14C SD"),
caption = "Per-layer soil observations from CiPEHR"
)
| Year | Layer | C stock (kg/m2) | C stock SD | Delta14C | Delta14C SD |
|---|---|---|---|---|---|
| 2009 | Surface organic | 1.9 | 0.5 | 151.0 | 56.5 |
| 2009 | Deep organic | 10.8 | 2.6 | -59.8 | 129.0 |
| 2009 | Mineral | 14.3 | 7.8 | -240.0 | 63.9 |
| 2022 | Surface organic | 1.9 | 0.5 | 136.2 | 62.7 |
| 2022 | Deep organic | 10.8 | 2.4 | -20.1 | 22.5 |
| 2022 | Mineral | 11.8 | 7.9 | -255.8 | 89.2 |
p_layer_c14 <- ggplot(dat_layers, aes(x = year, y = C14,
color = layer_label, shape = layer_label)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = C14 - C_14_sd, ymax = C14 + C_14_sd), width = 0.3) +
geom_line(linewidth = 0.8) +
labs(title = expression(paste(Delta^14, "C by Soil Layer")),
x = "Year", y = expression(paste(Delta^14, "C (per mil)")),
color = NULL, shape = NULL) +
scale_color_manual(values = c("coral", "steelblue", "forestgreen")) +
theme_minimal(base_size = 14) +
theme(legend.position = "top")
p_layer_cstock <- ggplot(dat_layers, aes(x = year, y = C_stock,
color = layer_label, shape = layer_label)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = C_stock - C_stock_sd, ymax = C_stock + C_stock_sd),
width = 0.3) +
geom_line(linewidth = 0.8) +
labs(title = "Carbon Stock by Soil Layer",
x = "Year", y = expression(paste("C stock (kg m"^-2, ")")),
color = NULL, shape = NULL) +
scale_color_manual(values = c("coral", "steelblue", "forestgreen")) +
theme_minimal(base_size = 14) +
theme(legend.position = "top")
p_layer_c14 + p_layer_cstock
Notice how the layers differ:
A 1-pool model can’t capture these differences. The 3-pool model treats each layer as a separate compartment, connected in series:
\[\text{Atmosphere} \xrightarrow{\text{input}} \text{Surface organic} \xrightarrow{a_{21}} \text{Deep organic} \xrightarrow{a_{32}} \text{Mineral}\]
The 3-pool series model has 6 parameters:
| Parameter | Description | Units |
|---|---|---|
k1 |
Decomposition rate, surface organic | yr\(^{-1}\) |
k2 |
Decomposition rate, deep organic | yr\(^{-1}\) |
k3 |
Decomposition rate, mineral | yr\(^{-1}\) |
input |
Carbon input rate to surface pool | gC/m\(^2\)/yr |
a21 |
Transfer rate, surface \(\to\) deep | yr\(^{-1}\) |
a32 |
Transfer rate, deep \(\to\) mineral | yr\(^{-1}\) |
The model follows the same structure as the 1-pool model – \(dC/dt = I + A \cdot C\) – but now \(C\) is a 3-element vector and \(A\) is a \(3 \times 3\) compartmental matrix:
\[A = \begin{pmatrix} -k_1 & 0 & 0 \\ a_{21} & -k_2 & 0 \\ 0 & a_{32} & -k_3 \end{pmatrix}\]
Each diagonal entry \(-k_i\) is the total loss rate from pool \(i\) (decomposition + transfer out). The off-diagonal entries \(a_{21}\) and \(a_{32}\) route carbon downward. This is a “series” topology – no carbon moves back up.
For this to be physically valid, the transfer rates must be smaller than the decomposition rates: \(a_{21} < k_1\) and \(a_{32} < k_2\). Otherwise, more carbon would leave a pool through transfer than through decomposition, which doesn’t make physical sense.
Let’s run the model at a few hand-picked parameter sets to build intuition.
dat_3pool <- datComb %>%
filter(year == 2009 | treatment == "Control") %>%
mutate(C_stock = C_stock * 1000,
C_stock_sd = C_stock_sd * 1000)
init_3pool <- dat_3pool %>% filter(year == 2009)
final_3pool <- dat_3pool %>% filter(year == 2022)
C0_3pool <- init_3pool$C_stock
F0_3pool <- init_3pool$C14
years_3pool <- seq(2009, 2022)
pool_labels <- c(so = "Surface organic", do = "Deep organic", min = "Mineral")
years_plot <- seq(2009, 2022, by = 0.5)
demo_params <- list(
list(ks = c(0.1, 0.01, 0.005), input = 100, a21 = 0.05, a32 = 0.005,
label = "Fast surface"),
list(ks = c(0.02, 0.01, 0.005), input = 50, a21 = 0.01, a32 = 0.002,
label = "Slow surface"),
list(ks = c(0.05, 0.02, 0.01), input = 80, a21 = 0.02, a32 = 0.01,
label = "High transfer")
)
demo_runs <- bind_rows(lapply(demo_params, function(dp) {
m <- tryCatch(
ThreepSeriesModel14(
t = years_plot,
ks = setNames(dp$ks, c("k1", "k2", "k3")),
C0 = C0_3pool, F0_Delta14C = F0_3pool,
In = dp$input, a21 = dp$a21, a32 = dp$a32,
inputFc = atmIn, pass = TRUE
),
error = function(e) NULL
)
if (is.null(m)) return(NULL)
data.frame(
year = rep(years_plot, 3),
pool = rep(pool_labels, each = length(years_plot)),
C14 = as.vector(getF14(m)),
C_stock = as.vector(getC(m)),
scenario = dp$label
)
}))
demo_runs$pool <- factor(demo_runs$pool, levels = pool_labels)
obs_plot <- dat_3pool %>%
mutate(pool = factor(pool_labels[layer], levels = pool_labels))
p_demo_c14 <- ggplot(demo_runs, aes(x = year, y = C14, color = scenario)) +
geom_line(linewidth = 1) +
geom_point(data = obs_plot, aes(x = year, y = C14),
inherit.aes = FALSE, size = 3) +
geom_errorbar(data = obs_plot,
aes(x = year, ymin = C14 - C_14_sd, ymax = C14 + C_14_sd),
inherit.aes = FALSE, width = 0.3) +
facet_wrap(~ pool, ncol = 1, scales = "free_y") +
labs(title = expression(paste("Forward model demo: ", Delta^14, "C")),
x = "Year", y = expression(paste(Delta^14, "C (per mil)")),
color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top")
p_demo_cstock <- ggplot(demo_runs, aes(x = year, y = C_stock / 1000, color = scenario)) +
geom_line(linewidth = 1) +
geom_point(data = obs_plot, aes(x = year, y = C_stock / 1000),
inherit.aes = FALSE, size = 3) +
geom_errorbar(data = obs_plot,
aes(x = year, ymin = (C_stock - C_stock_sd) / 1000,
ymax = (C_stock + C_stock_sd) / 1000),
inherit.aes = FALSE, width = 0.3) +
facet_wrap(~ pool, ncol = 1, scales = "free_y") +
labs(title = "Forward model demo: Carbon stock",
x = "Year", y = expression(paste("C stock (kg m"^-2, ")")),
color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top")
p_demo_c14 + p_demo_cstock
None of these hand-picked parameter sets perfectly match the observations in all three layers simultaneously. That’s the challenge – and exactly why we need MCMC to explore the 6-dimensional parameter space.
The likelihood has the same Gaussian structure as the 1-pool version, but now we compare predictions to observations in each of the 3 pools separately. That gives us 6 data points: 3 carbon stocks and 3 \(\Delta^{14}\)C values at 2022.
One important difference: we multiply the reported measurement uncertainty by a factor of 0.4. This is an empirical scaling factor from the original analysis that accounts for model structural error – the fact that the model is a simplified representation of reality. Without this scaling, the likelihood would be too forgiving and the posteriors too wide.
likelihood_3pool <- function(pars) {
k1 <- pars[1]; k2 <- pars[2]; k3 <- pars[3]
input <- pars[4]; a21 <- pars[5]; a32 <- pars[6]
model <- tryCatch(
ThreepSeriesModel14(
t = years_3pool,
ks = c(k1 = k1, k2 = k2, k3 = k3),
C0 = C0_3pool,
F0_Delta14C = F0_3pool,
In = input,
a21 = a21,
a32 = a32,
inputFc = atmIn,
pass = TRUE
),
error = function(e) return(NULL)
)
if (is.null(model)) return(-1e10)
idx <- length(years_3pool)
pred_C14 <- getF14(model)[idx, ]
pred_Cstock <- getC(model)[idx, ]
ll_C14 <- sum(dnorm(pred_C14 - final_3pool$C14, 0,
final_3pool$C_14_sd * 0.4, log = TRUE))
ll_Cstock <- sum(dnorm(pred_Cstock - final_3pool$C_stock, 0,
final_3pool$C_stock_sd * 0.4, log = TRUE))
return(sum(ll_C14, ll_Cstock, na.rm = TRUE))
}
cat("Test likelihood at example parameters:",
likelihood_3pool(c(0.1, 0.01, 0.01, 200, 0.1, 0.01)))
## Test likelihood at example parameters: -42.98645
Unlike the 1-pool model, we can’t easily visualize a 6-dimensional likelihood surface. We have to trust the MCMC sampler to explore it for us. This is where good priors and sufficient iterations become critical.
As we discussed in the 1-pool tutorial, the prior encodes what we believe before seeing the data. With 6 parameters, prior choice matters even more – too wide and the MCMC wastes time in impossible regions; too narrow and you might miss the answer.
We have two broad options:
We’ll start with uniform priors and then compare to truncated normal priors in Section 9.
The key physical constraints for both prior types:
par_names <- c("k1", "k2", "k3", "input", "a21", "a32")
prior_lower <- c(k1 = 0.002, k2 = 0.0006, k3 = 0.0002,
input = 10, a21 = 0.0025, a32 = 0.0001)
prior_upper <- c(k1 = 0.2, k2 = 0.06, k3 = 0.06,
input = 205, a21 = 0.5, a32 = 0.1)
prior_range_df <- data.frame(
Parameter = par_names,
Lower = prior_lower,
Upper = prior_upper,
Interpretation = c(
"Turnover 5-500 yr",
"Turnover 17-1667 yr",
"Turnover 17-5000 yr",
"Litter input to surface pool",
"Surface to deep transfer",
"Deep to mineral transfer"
)
)
knitr::kable(prior_range_df, digits = 4,
caption = "Uniform prior ranges for the 3-pool model")
| Parameter | Lower | Upper | Interpretation | |
|---|---|---|---|---|
| k1 | k1 | 0.0020 | 0.20 | Turnover 5-500 yr |
| k2 | k2 | 0.0006 | 0.06 | Turnover 17-1667 yr |
| k3 | k3 | 0.0002 | 0.06 | Turnover 17-5000 yr |
| input | input | 10.0000 | 205.00 | Litter input to surface pool |
| a21 | a21 | 0.0025 | 0.50 | Surface to deep transfer |
| a32 | a32 | 0.0001 | 0.10 | Deep to mineral transfer |
Suppose a literature review suggests approximate values for each parameter. We can encode this as truncated normal priors – Gaussian bell curves centered on the literature values, clipped to the same physical bounds as the uniform priors.
# Literature-informed centers and uncertainties
prior_means <- c(k1 = 0.05, k2 = 0.015, k3 = 0.005,
input = 80, a21 = 0.03, a32 = 0.01)
prior_sds <- c(k1 = 0.04, k2 = 0.015, k3 = 0.01,
input = 60, a21 = 0.05, a32 = 0.02)
prior_gaussian <- createTruncatedNormalPrior(
mean = prior_means,
sd = prior_sds,
lower = prior_lower,
upper = prior_upper
)
gaussian_range_df <- data.frame(
Parameter = par_names,
Mean = prior_means,
SD = prior_sds,
Lower = prior_lower,
Upper = prior_upper
)
knitr::kable(gaussian_range_df, digits = 4,
caption = "Truncated normal prior parameters (same bounds as uniform)")
| Parameter | Mean | SD | Lower | Upper | |
|---|---|---|---|---|---|
| k1 | k1 | 0.050 | 0.040 | 0.0020 | 0.20 |
| k2 | k2 | 0.015 | 0.015 | 0.0006 | 0.06 |
| k3 | k3 | 0.005 | 0.010 | 0.0002 | 0.06 |
| input | input | 80.000 | 60.000 | 10.0000 | 205.00 |
| a21 | a21 | 0.030 | 0.050 | 0.0025 | 0.50 |
| a32 | a32 | 0.010 | 0.020 | 0.0001 | 0.10 |
# Uniform prior samples
prior_samples_unif <- as.data.frame(mapply(
runif, n = 10000, min = prior_lower, max = prior_upper
))
# Truncated normal prior samples
prior_samples_gauss <- as.data.frame(mapply(function(mu, s, lo, hi) {
# Rejection sampling from truncated normal
out <- numeric(10000)
filled <- 0
while (filled < 10000) {
candidates <- rnorm(20000, mu, s)
valid <- candidates[candidates >= lo & candidates <= hi]
n_take <- min(length(valid), 10000 - filled)
out[(filled + 1):(filled + n_take)] <- valid[1:n_take]
filled <- filled + n_take
}
out
}, mu = prior_means, s = prior_sds, lo = prior_lower, hi = prior_upper))
colnames(prior_samples_gauss) <- par_names
prior_long_unif <- prior_samples_unif %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
mutate(parameter = factor(parameter, levels = par_names),
prior_type = "Uniform")
prior_long_gauss <- prior_samples_gauss %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
mutate(parameter = factor(parameter, levels = par_names),
prior_type = "Gaussian")
prior_long <- bind_rows(prior_long_unif, prior_long_gauss)
ggplot(prior_long, aes(x = value, fill = prior_type)) +
geom_density(alpha = 0.4, linewidth = 0.6) +
facet_wrap(~ parameter, scales = "free", ncol = 3) +
scale_fill_manual(values = c("Gaussian" = "coral", "Uniform" = "skyblue"),
name = "Prior type") +
labs(title = "Uniform vs. Gaussian (Truncated Normal) Priors",
subtitle = "Same physical bounds, but Gaussian concentrates probability near the center",
x = "Parameter value", y = "Density") +
theme_minimal(base_size = 12) +
theme(legend.position = "top")
The Gaussian priors have the same physical bounds as the uniform priors, but they concentrate probability near the literature-informed centers. For well-constrained parameters (like \(k_1\)), this may not matter much – the data will drive the posterior regardless. For weakly constrained parameters (like \(k_3\) or \(a_{32}\)), the prior shape can make a real difference.
With 6 parameters, we need more iterations than the 1-pool model. We use 20,000 iterations with DREAMzs. For a teaching example this is sufficient; a publication-quality analysis might use 100,000+.
result <- run_and_summarize(
likelihood_fn = likelihood_3pool,
prior_lower = prior_lower,
prior_upper = prior_upper,
par_names = par_names,
iterations = 20000,
burn_in = 5000,
thin = 10
)
cat("Retained", nrow(result$posterior_df), "posterior samples")
## Retained 501 posterior samples
With 6 parameters, we need a bigger panel. Look for the same “hairy caterpillar” pattern in each – well-mixed chains bouncing around a stable region.
trace_vars <- c(par_names, "log_posterior")
trace_labels <- c(par_names, "Log-posterior")
chain_long <- result$chain_df %>%
pivot_longer(cols = all_of(c(par_names, "log_posterior")),
names_to = "parameter", values_to = "value") %>%
mutate(parameter = factor(parameter, levels = trace_vars, labels = trace_labels))
ggplot(chain_long, aes(x = iteration, y = value, color = chain)) +
geom_line(alpha = 0.4, linewidth = 0.3) +
facet_wrap(~ parameter, scales = "free_y", ncol = 3) +
labs(title = "MCMC Trace Plots (3-Pool Model)",
subtitle = "Log-posterior (bottom right) should plateau when chains converge",
x = "Iteration", y = "Value") +
theme_minimal(base_size = 11) +
theme(legend.position = "top")
The log-posterior panel is the most important diagnostic: it shows the objective value being maximized. Look for a clear plateau – the point where the log-posterior stops climbing and starts bouncing around a stable level. This is where the chains have found the high-probability region, and everything before it is burn-in.
Some parameters (like k3 or a32 for the
slow mineral pool) may show worse mixing than others. This is typical –
slow pools have weak 14C signals, so the data provide less constraint
and the chains wander more.
We use a burn-in of 5,000 (vs. 2,000 for the 1-pool model) because the chains need more time to settle in a 6D space. Thinning of 10 reduces autocorrelation.
ggplot(chain_long, aes(x = iteration, y = value, color = chain)) +
geom_line(alpha = 0.4, linewidth = 0.3) +
geom_vline(xintercept = result$burn_in / 3, color = "black",
linewidth = 1, linetype = "dotted") +
annotate("text", x = result$burn_in / 3, y = Inf,
label = " Burn-in cutoff", hjust = 0, vjust = 1.5, size = 3.5) +
facet_wrap(~ parameter, scales = "free_y", ncol = 3) +
labs(title = "Burn-in Removal",
subtitle = "Everything left of the dotted line is discarded. Check that log-posterior has plateaued by the cutoff.",
x = "Iteration", y = "Value") +
theme_minimal(base_size = 11) +
theme(legend.position = "top")
The key question: how much did the data narrow down our estimates for each parameter?
posterior_long <- result$posterior_df %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
mutate(parameter = factor(parameter, levels = par_names))
ggplot() +
geom_histogram(data = prior_long_unif,
aes(x = value, y = after_stat(density)),
bins = 40, fill = "skyblue", alpha = 0.4, color = "white") +
geom_density(data = posterior_long, aes(x = value),
fill = "coral", alpha = 0.5, linewidth = 0.8) +
geom_vline(data = data.frame(
parameter = factor(par_names, levels = par_names),
xval = result$map),
aes(xintercept = xval), linetype = "dashed", color = "darkblue") +
facet_wrap(~ parameter, scales = "free", ncol = 3) +
labs(title = "Prior (blue) vs. Posterior (coral)",
subtitle = "Dashed blue = MAP estimate. Narrow posteriors = well-constrained parameters.",
x = "Parameter value", y = "Density") +
theme_minimal(base_size = 12)
# Compute uncertainty reduction
summary_with_reduction <- result$summary %>%
mutate(
prior_range = prior_upper[parameter] - prior_lower[parameter],
post_range = upper_90 - lower_90,
uncert_reduction = round(100 * (1 - post_range / prior_range), 1)
)
cat("=== 3-Pool Posterior Summary ===\n")
## === 3-Pool Posterior Summary ===
cat(sprintf("%-8s %8s %10s %10s %10s %s\n",
"Param", "MAP", "Median", "90% lo", "90% hi", "Reduction"))
## Param MAP Median 90% lo 90% hi Reduction
for (i in seq_len(nrow(summary_with_reduction))) {
r <- summary_with_reduction[i, ]
cat(sprintf("%-8s %8.4f %10.4f %10.4f %10.4f %5.1f%%\n",
r$parameter, result$map[r$parameter], r$median,
r$lower_90, r$upper_90, r$uncert_reduction))
}
## a21 0.0902 0.1057 0.0695 0.1573 82.4%
## a32 0.0022 0.0119 0.0015 0.0317 69.8%
## input 30.4176 46.9212 16.7631 106.4977 54.0%
## k1 0.0221 0.0323 0.0081 0.0635 72.0%
## k2 0.0165 0.0188 0.0060 0.0348 51.4%
## k3 0.0123 0.0256 0.0040 0.0518 20.1%
The turnover time \(\tau = 1/k\) tells us how long carbon resides in each pool on average. This is often more intuitive than the decomposition rate.
k_params <- c("k1", "k2", "k3")
turnover_df <- result$posterior_df[, k_params] %>%
mutate(across(everything(), ~ 1 / .x)) %>%
rename(tau1 = k1, tau2 = k2, tau3 = k3) %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
mutate(pool = factor(parameter,
levels = c("tau1", "tau2", "tau3"),
labels = c("Surface organic", "Deep organic", "Mineral")))
turnover_summary <- turnover_df %>%
group_by(pool) %>%
summarise(
median = median(value),
lower_90 = quantile(value, 0.05),
upper_90 = quantile(value, 0.95),
.groups = "drop"
)
knitr::kable(turnover_summary, digits = 1,
col.names = c("Pool", "Median (yr)", "90% lo", "90% hi"),
caption = "Posterior turnover times (1/k)")
| Pool | Median (yr) | 90% lo | 90% hi |
|---|---|---|---|
| Surface organic | 30.9 | 15.7 | 123.2 |
| Deep organic | 53.2 | 28.7 | 167.4 |
| Mineral | 39.1 | 19.3 | 248.5 |
In multi-parameter models, parameters often trade off against each other. For example, a higher decomposition rate might be compensated by higher carbon input. These correlations are a key feature of the posterior that single-parameter summaries miss.
cor_mat <- cor(result$posterior_df)
cat("--- Posterior Correlation Matrix ---\n")
## --- Posterior Correlation Matrix ---
print(round(cor_mat, 2))
## k1 k2 k3 input a21 a32
## k1 1.00 0.14 0.09 0.84 0.33 -0.29
## k2 0.14 1.00 0.01 0.22 0.42 0.03
## k3 0.09 0.01 1.00 0.09 0.05 0.24
## input 0.84 0.22 0.09 1.00 0.27 -0.23
## a21 0.33 0.42 0.05 0.27 1.00 -0.09
## a32 -0.29 0.03 0.24 -0.23 -0.09 1.00
n_pairs <- min(500, nrow(result$posterior_df))
pairs(result$posterior_df[sample(nrow(result$posterior_df), n_pairs), ],
pch = 16, cex = 0.4, col = adjustcolor("steelblue", 0.3),
main = "Posterior Pairwise Correlations (3-Pool Model)")
Look for elliptical clouds tilted at an angle – these indicate correlated parameters. Strong correlations mean the data can’t independently constrain those parameters; only certain combinations are well determined.
We run the forward model at the MAP estimate (best fit) and at 200 random posterior draws (uncertainty envelope), showing predictions for each soil layer separately.
years_plot <- seq(2009, 2022, by = 0.5)
# MAP prediction
map_model <- ThreepSeriesModel14(
t = years_plot,
ks = c(k1 = result$map["k1"], k2 = result$map["k2"], k3 = result$map["k3"]),
C0 = C0_3pool, F0_Delta14C = F0_3pool,
In = result$map["input"], a21 = result$map["a21"], a32 = result$map["a32"],
inputFc = atmIn, pass = TRUE
)
map_pred <- data.frame(
year = rep(years_plot, 3),
pool = rep(pool_labels, each = length(years_plot)),
C14 = as.vector(getF14(map_model)),
C_stock = as.vector(getC(map_model))
)
# Posterior ensemble
n_draws <- min(200, nrow(result$posterior_df))
draw_idx <- sample(nrow(result$posterior_df), n_draws)
ensemble <- bind_rows(lapply(draw_idx, function(i) {
p <- as.numeric(result$posterior_df[i, ])
m <- tryCatch(
ThreepSeriesModel14(
t = years_plot,
ks = c(k1 = p[1], k2 = p[2], k3 = p[3]),
C0 = C0_3pool, F0_Delta14C = F0_3pool,
In = p[4], a21 = p[5], a32 = p[6],
inputFc = atmIn, pass = TRUE
),
error = function(e) NULL
)
if (is.null(m)) return(NULL)
data.frame(
year = rep(years_plot, 3),
pool = rep(pool_labels, each = length(years_plot)),
C14 = as.vector(getF14(m)),
C_stock = as.vector(getC(m))
)
}))
envelope <- ensemble %>%
group_by(year, pool) %>%
summarise(
C14_lo = quantile(C14, 0.05), C14_hi = quantile(C14, 0.95),
Cstock_lo = quantile(C_stock, 0.05), Cstock_hi = quantile(C_stock, 0.95),
.groups = "drop"
)
obs_plot <- dat_3pool %>%
mutate(pool = factor(pool_labels[layer], levels = pool_labels))
map_pred$pool <- factor(map_pred$pool, levels = pool_labels)
envelope$pool <- factor(envelope$pool, levels = pool_labels)
# Delta14C by layer
p_c14 <- ggplot() +
geom_ribbon(data = envelope,
aes(x = year, ymin = C14_lo, ymax = C14_hi),
fill = "steelblue", alpha = 0.3) +
geom_line(data = map_pred, aes(x = year, y = C14),
color = "steelblue", linewidth = 1) +
geom_point(data = obs_plot, aes(x = year, y = C14), size = 3) +
geom_errorbar(data = obs_plot,
aes(x = year, ymin = C14 - C_14_sd, ymax = C14 + C_14_sd),
width = 0.3) +
facet_wrap(~ pool, ncol = 1, scales = "free_y") +
labs(title = expression(paste(Delta^14, "C by Soil Layer")),
subtitle = "MAP prediction + 90% posterior envelope vs. observations",
x = "Year", y = expression(paste(Delta^14, "C (per mil)"))) +
theme_minimal(base_size = 12)
# C stock by layer
p_cstock <- ggplot() +
geom_ribbon(data = envelope,
aes(x = year, ymin = Cstock_lo / 1000, ymax = Cstock_hi / 1000),
fill = "forestgreen", alpha = 0.3) +
geom_line(data = map_pred, aes(x = year, y = C_stock / 1000),
color = "forestgreen", linewidth = 1) +
geom_point(data = obs_plot, aes(x = year, y = C_stock / 1000), size = 3) +
geom_errorbar(data = obs_plot,
aes(x = year, ymin = (C_stock - C_stock_sd) / 1000,
ymax = (C_stock + C_stock_sd) / 1000),
width = 0.3) +
facet_wrap(~ pool, ncol = 1, scales = "free_y") +
labs(title = "Carbon Stock by Soil Layer",
subtitle = "MAP prediction + 90% posterior envelope vs. observations",
x = "Year", y = expression(paste("C stock (kg m"^-2, ")"))) +
theme_minimal(base_size = 12)
p_c14 + p_cstock
Does prior shape matter for the 3-pool model? Let’s re-run the MCMC with the truncated normal priors from Section 5 and compare the posteriors.
result_gauss <- run_and_summarize(
likelihood_fn = likelihood_3pool,
prior = prior_gaussian,
par_names = par_names,
iterations = 20000,
burn_in = 5000,
thin = 10
)
cat("Retained", nrow(result_gauss$posterior_df), "posterior samples (Gaussian prior)")
## Retained 501 posterior samples (Gaussian prior)
posterior_long_unif <- result$posterior_df %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
mutate(parameter = factor(parameter, levels = par_names),
prior_type = "Uniform")
posterior_long_gauss <- result_gauss$posterior_df %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value") %>%
mutate(parameter = factor(parameter, levels = par_names),
prior_type = "Gaussian")
posteriors_compare <- bind_rows(posterior_long_unif, posterior_long_gauss)
ggplot(posteriors_compare, aes(x = value, fill = prior_type)) +
geom_density(alpha = 0.4, linewidth = 0.6) +
facet_wrap(~ parameter, scales = "free", ncol = 3) +
scale_fill_manual(values = c("Gaussian" = "coral", "Uniform" = "steelblue"),
name = "Prior type") +
labs(title = "Posterior Comparison: Uniform vs. Gaussian Priors",
subtitle = "Same likelihood and data, different prior shapes",
x = "Parameter value", y = "Density") +
theme_minimal(base_size = 12) +
theme(legend.position = "top")
summary_gauss <- result_gauss$summary %>%
mutate(
prior_range = prior_upper[parameter] - prior_lower[parameter],
post_range = upper_90 - lower_90,
uncert_reduction = round(100 * (1 - post_range / prior_range), 1)
)
cat("=== Uniform Prior: Posterior Summary ===\n")
## === Uniform Prior: Posterior Summary ===
cat(sprintf("%-8s %10s %10s %10s\n", "Param", "Median", "90% lo", "90% hi"))
## Param Median 90% lo 90% hi
for (i in seq_len(nrow(summary_with_reduction))) {
r <- summary_with_reduction[i, ]
cat(sprintf("%-8s %10.4f %10.4f %10.4f\n",
r$parameter, r$median, r$lower_90, r$upper_90))
}
## a21 0.1057 0.0695 0.1573
## a32 0.0119 0.0015 0.0317
## input 46.9212 16.7631 106.4977
## k1 0.0323 0.0081 0.0635
## k2 0.0188 0.0060 0.0348
## k3 0.0256 0.0040 0.0518
cat("\n=== Gaussian Prior: Posterior Summary ===\n")
##
## === Gaussian Prior: Posterior Summary ===
cat(sprintf("%-8s %10s %10s %10s\n", "Param", "Median", "90% lo", "90% hi"))
## Param Median 90% lo 90% hi
for (i in seq_len(nrow(summary_gauss))) {
r <- summary_gauss[i, ]
cat(sprintf("%-8s %10.4f %10.4f %10.4f\n",
r$parameter, r$median, r$lower_90, r$upper_90))
}
## a21 0.0907 0.0566 0.1285
## a32 0.0078 0.0009 0.0261
## input 50.1181 18.6064 106.6724
## k1 0.0307 0.0123 0.0553
## k2 0.0180 0.0058 0.0308
## k3 0.0115 0.0018 0.0236
What to notice:
The practical takeaway: when data are limited (as they often are in soil science), the shape of your prior is part of your scientific argument. Gaussian priors informed by literature values are a principled way to incorporate existing knowledge. But always run a sensitivity comparison like this one to verify that your conclusions aren’t driven by the prior alone.
How does the 3-pool model compare to the 1-pool model? Let’s run a quick 1-pool calibration for comparison and look at bulk predictions from both models.
# Aggregate to bulk observations (same as 1-pool tutorial)
dat_1pool <- datComb %>%
filter(year == 2009 | treatment == "Control") %>%
mutate(C_stock_g = C_stock * 1000,
C_stock_sd_g = C_stock_sd * 1000) %>%
group_by(year) %>%
summarise(
C14 = weighted.mean(C14, C_stock_g),
C_14_sd = sqrt(sum((C_stock_g * C_14_sd)^2)) / sum(C_stock_g),
C_stock = sum(C_stock_g),
C_stock_sd = sqrt(sum(C_stock_sd_g^2)),
.groups = "drop"
)
C0_1pool <- dat_1pool$C_stock[dat_1pool$year == 2009]
F0_1pool <- dat_1pool$C14[dat_1pool$year == 2009]
years_1pool <- seq(2009, 2022)
obs_2022_1pool <- dat_1pool %>% filter(year == 2022)
likelihood_1pool <- function(pars) {
k <- pars[1]
input <- pars[2]
model <- tryCatch(
OnepModel14(t = years_1pool, k = k, C0 = C0_1pool,
F0_Delta14C = F0_1pool, In = input,
inputFc = atmIn, pass = TRUE),
error = function(e) return(NULL)
)
if (is.null(model)) return(-1e10)
pred_C14 <- getF14(model)[length(years_1pool), ]
pred_Cstock <- getC(model)[length(years_1pool), ]
ll_C14 <- dnorm(pred_C14 - obs_2022_1pool$C14, 0,
obs_2022_1pool$C_14_sd, log = TRUE)
ll_Cstock <- dnorm(pred_Cstock - obs_2022_1pool$C_stock, 0,
obs_2022_1pool$C_stock_sd, log = TRUE)
return(sum(ll_C14, ll_Cstock, na.rm = TRUE))
}
result_1pool <- run_and_summarize(
likelihood_fn = likelihood_1pool,
prior_lower = c(k = 0.001, input = 50),
prior_upper = c(k = 0.1, input = 1000),
par_names = c("k", "input"),
iterations = 10000,
burn_in = 2000,
thin = 5
)
years_plot_cmp <- seq(2009, 2022, by = 0.5)
# 1-pool MAP prediction
map_1pool <- OnepModel14(
t = years_plot_cmp, k = result_1pool$map["k"], C0 = C0_1pool,
F0_Delta14C = F0_1pool, In = result_1pool$map["input"],
inputFc = atmIn, pass = TRUE
)
pred_1pool <- data.frame(
year = years_plot_cmp,
C14 = as.numeric(getF14(map_1pool)),
C_stock = as.numeric(getC(map_1pool)),
model = "1-pool"
)
# 3-pool MAP prediction, aggregated to bulk
map_3pool_model <- ThreepSeriesModel14(
t = years_plot_cmp,
ks = c(k1 = result$map["k1"], k2 = result$map["k2"], k3 = result$map["k3"]),
C0 = C0_3pool, F0_Delta14C = F0_3pool,
In = result$map["input"], a21 = result$map["a21"], a32 = result$map["a32"],
inputFc = atmIn, pass = TRUE
)
C_3pool <- getC(map_3pool_model)
C14_3pool <- getF14(map_3pool_model)
# Bulk = sum of C stocks; bulk 14C = stock-weighted mean
bulk_C <- rowSums(C_3pool)
bulk_C14 <- rowSums(C_3pool * C14_3pool) / bulk_C
pred_3pool <- data.frame(
year = years_plot_cmp,
C14 = as.numeric(bulk_C14),
C_stock = as.numeric(bulk_C),
model = "3-pool (bulk)"
)
pred_compare <- bind_rows(pred_1pool, pred_3pool)
p_cmp_c14 <- ggplot(pred_compare, aes(x = year, y = C14, color = model)) +
geom_line(linewidth = 1.2) +
geom_point(data = dat_1pool, aes(x = year, y = C14),
inherit.aes = FALSE, size = 4) +
geom_errorbar(data = dat_1pool,
aes(x = year, ymin = C14 - C_14_sd, ymax = C14 + C_14_sd),
inherit.aes = FALSE, width = 0.3) +
scale_color_manual(values = c("steelblue", "coral")) +
labs(title = expression(paste("Bulk ", Delta^14, "C: 1-Pool vs. 3-Pool")),
x = "Year", y = expression(paste(Delta^14, "C (per mil)")),
color = NULL) +
theme_minimal(base_size = 14) +
theme(legend.position = "top")
p_cmp_cstock <- ggplot(pred_compare, aes(x = year, y = C_stock / 1000, color = model)) +
geom_line(linewidth = 1.2) +
geom_point(data = dat_1pool, aes(x = year, y = C_stock / 1000),
inherit.aes = FALSE, size = 4) +
geom_errorbar(data = dat_1pool,
aes(x = year, ymin = (C_stock - C_stock_sd) / 1000,
ymax = (C_stock + C_stock_sd) / 1000),
inherit.aes = FALSE, width = 0.3) +
scale_color_manual(values = c("steelblue", "coral")) +
labs(title = "Bulk Carbon Stock: 1-Pool vs. 3-Pool",
x = "Year", y = expression(paste("C stock (kg m"^-2, ")")),
color = NULL) +
theme_minimal(base_size = 14) +
theme(legend.position = "top")
p_cmp_c14 / p_cmp_cstock
What does the 3-pool model tell us that the 1-pool model can’t?
Layer-specific turnover times: The 3-pool model estimates that surface organic carbon cycles in ~45 years, while mineral carbon takes ~81 years. The 1-pool model gives a single bulk estimate that averages over these very different timescales.
Carbon transfer through the profile: The transfer rates \(a_{21}\) and \(a_{32}\) quantify how fast carbon moves downward from surface to deep to mineral soil. This is invisible to the 1-pool model.
Layer-specific predictions: The 3-pool model predicts how each individual layer will evolve, not just the bulk. This matters for understanding vulnerability – surface organic carbon may respond very differently to warming than deep mineral carbon.
The cost of complexity: More parameters mean wider posteriors (less certainty per parameter), stronger correlations between parameters, and harder convergence. The 6-parameter posterior is much harder for MCMC to explore than the 2-parameter posterior. This is the classic bias-variance tradeoff in modeling – more complex models can capture more reality, but they’re harder to constrain from limited data.
Multi-pool models capture real structure that single-pool models miss. The three soil layers have genuinely different turnover times, and a 3-pool model can estimate these separately.
More parameters = harder inference. With 6 parameters and 6 data points, parameter correlations become important. Always check the posterior correlation matrix – it tells you which parameters the data can and can’t independently constrain.
Prior shape matters for weakly constrained parameters. We saw in Section 9 that switching from uniform to Gaussian priors can noticeably affect posteriors for parameters with weak data constraint. When data are limited, the prior is part of your scientific argument – choose it deliberately and always run a sensitivity analysis.
Trace plot diagnosis is even more critical with more parameters. Some parameters (especially for slow pools) may mix poorly because the data provide weak constraints. This is a feature, not a bug – it tells you where the data are uninformative.
The uncertainty scaling factor (0.4) on the likelihood matters a lot. It accounts for the gap between measurement uncertainty and model structural error. Experiment with it to understand its effect.
Model comparison is essential. Comparing 1-pool and 3-pool predictions side by side (Section 10) reveals when the extra complexity is justified and when a simpler model suffices.
Try feedback topology: Instead of a pure series model (carbon only flows downward), allow some upward transfer. How does this change the posteriors?
Try the warming treatment: The
datComb table also contains a “Warming” treatment. What
happens to decomposition rates under experimental warming?
Try log-normal priors for k: Decomposition rates
span orders of magnitude, making log-normal priors a natural choice. Use
createPrior() with a custom density and sampler
function.
Increase iterations: Try 100,000 iterations. Do the posteriors change? If so, 20,000 wasn’t enough for full convergence.
sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Phoenix
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.3.2 tidyr_1.3.1 dplyr_1.1.4
## [4] ggplot2_4.0.0 BayesianTools_0.1.8 SoilR_1.2.107
## [7] deSolve_1.40
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.52 bslib_0.7.0
## [4] sets_1.0-25 lattice_0.22-6 vctrs_0.6.5
## [7] tools_4.4.0 Rdpack_2.6.6 generics_0.1.3
## [10] parallel_4.4.0 tibble_3.2.1 fansi_1.0.6
## [13] highr_0.10 pkgconfig_2.0.3 Matrix_1.7-0
## [16] DHARMa_0.4.7 RColorBrewer_1.1-3 S7_0.2.0
## [19] assertthat_0.2.1 lifecycle_1.0.4 compiler_4.4.0
## [22] farver_2.1.2 stringr_1.5.1 Brobdingnag_1.2-9
## [25] htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.8
## [28] pillar_1.9.0 nloptr_2.2.1 jquerylib_0.1.4
## [31] MASS_7.3-60.2 cachem_1.1.0 reformulas_0.4.4
## [34] bridgesampling_1.2-1 boot_1.3-30 nlme_3.1-164
## [37] tidyselect_1.2.1 digest_0.6.35 mvtnorm_1.2-4
## [40] stringi_1.8.4 purrr_1.1.0 labeling_0.4.3
## [43] splines_4.4.0 fastmap_1.2.0 grid_4.4.0
## [46] expm_1.0-0 cli_3.6.2 magrittr_2.0.3
## [49] survival_3.5-8 utf8_1.2.4 withr_3.0.0
## [52] scales_1.4.0 rmarkdown_2.29 igraph_2.1.4
## [55] lme4_1.1-38 msm_1.8.2 coda_0.19-4.1
## [58] evaluate_0.23 knitr_1.46 rbibutils_2.4.1
## [61] rlang_1.1.3 Rcpp_1.0.12 glue_1.7.0
## [64] minqa_1.2.8 jsonlite_1.8.8 R6_2.5.1